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ABSTRACT 


A recently developed finitc-elcmcnt-based equivalent linearization approach 
for the analysis of random vibrations of geometrically nonlinear multiple 
degree-of-freedom structures is validated. The validation is based on 
comparisons with results from a finite element based numerical simulation 
analysis using a numerical integration technique in physical coordinates. In 
particular, results for the case of a clamped-clamped beam are considered for 
an extensive load range to establish the limits of validity of the equivalent 
linearization approach. 

INTRODUCTION 

Current efforts to extend the perfonnance and flight envelope of high-speed 
aerospace vehicles have resulted in structures which may respond to the 
imposed loads in a geometrically nonlinear (large deflection) random fashion. 
This type of behavior can significantly degrade the structural fatigue life. 
Linear prediction techniques currently used in the design process are grossly 
conservative and provide little understanding of the nonlinear behavior. 
Without practical design tools capable of capturing the important dynamics, 
further improvements in vehicle performance and system design will be 
hampered. 

Methods currently used to predict geometrically nonlinear random response 
include perturbation, Fokker-Plank-Kolmogorov (F-P-K), numerical 
simulation and stochastic linearization techniques. All have various 
limitations. Perturbation techniques are limited to weak geometric 
nonlinearities. The F-P-K approach [1,2] yields exact solutions, but can only 
be applied to simple mechanical systems. Numerical simulation techniques 
using numerical integration provide time histories of the response from which 
statistics of the random response may be calculated. This, however, comes at 
a high computational expense due to the long time records or high number of 
ensemble averages required to get quality random response statistics. 
Statistical linearization methods (e.g. equivalent linearization (EL), see [2-6]) 
have seen the most broad application because of their ability to accurately 
capture the response statistics over a wide range of response levels while 
maintaining relatively light computational burden. Implementations of 
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statistical linearization methods have been primarily limited to special purpose 
“in-house” computer codes. The first known implementation in a general- 
purpose code was developed in [7] using MSC/NASTRAN. In this study, an 
alternative implementation of EL based on the methods developed in [8] and 
[9] will be used. This implementation was previously validated using the F-P- 
K method for a few special cases including a Duffing oscillator and beam 
structure under a convenient, but non-physical loading condition. In the 
present study, results from the EL analysis will be compared with those from a 
finite element based numerical simulation analysis for the case of a clamped- 
clamped beam under random inertial loading. By studying a wide range of 
load levels, the range of applicability is established. 

EQUIVALENT LINEARIZATION APPROACH 

The equations of motion of a multiple degree-of-freedom, viscously damped 
geometrically nonlinear system can be written in the form: 

MX(t) + CX(t ) + KX(t) + r(X(t)) = F{t) (1) 


where M, C, K are the mass, damping, and stiffness matrices, X is the 
displacement response vector and F is the force excitation vector, 
respectively. The nonlinear stiffness tenn T(X) is a vector function which 
gen erally includes 2 nd and 3 rd order terms in X . An approximate solution to 
(1) can be achieved by formation of an equivalent linear system: 

MX(t) + CX(t) + (K + K e )X(t) = F(t ) (2) 


where K e is the equivalent linear stiffness matrix. 


The traditional (force error minimization) method of EL seeks to minimize the 
difference between the nonlinear force and the product of the equivalent linear 
stiffness and displacement response vector. Since the error is a random 
function of time, the required condition is that the expectation of the mean 
square error be a minimum, i.e.: 

error — E[(T(X) - K e Xf( T(X) - KX)\ -> min (3) 


where £[...] represents the expectation operator. Equation [(3)] will be satisfied 
if 


3 {error) _ 
~dK~~ 


ij = 1,2,..., A 


In this study, consideration is limited to the case of Gaussian, zero-mean 
excitation and response to simplify the solution. With these assumptions and 
omitting intermediate derivations, the final form for the equivalent linear 
stiffness matrix becomes (see for example [3] and [4]): 


K. 



( 4 ) 
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An alternative approach based on potential (strain) energy error minimization 
was proposed in [5] and [6], where mostly single degree-of- freedom systems 
were considered. A generalization for multiple degree-of-freedom systems is 
developed in [8] and [9]. For the sake of brevity, the present paper will present 
the formulation and results from the force error minimization only. 

Applying the modal coordinate transformation 

X - Oq 

to equation |(2)| yields a set of coupled modal equations with reduced degrees 
of freedom, where <f> is generally a subset (L< N) of the linear eigenvectors, 
and q are the modal coordinates. This coupled set is expressed as: 

UVl + [2 C r CO r ]q + [[©2] + [k^q = f (5) 

where co r are the undamped natural frequencies, [I] is the unit matrix, [2 £ r (O r \ 
is the diagonal modal damping matrix, [ftr] is the diagonal modal stiffness 
matrix, and [fcj is the fully populated equivalent stiffness matrix given by: 

<6> 

The nonlinear terms may be represented in the following form: 

L L 

7 r (qi , • • • , q L ) = X + X br j»<iflk<h ( 7 ) 

j,k=j j,k=j,l=k 

where a jk and b' jU are nonlinear stiffness coefficients with j = 1,2,..., L, k = j, 

j+l,...,L, and l = k, k+l, ...,L. This form of the nonlinear terms facilitates the 
solution of equations (5) when the forces and displacements are random 
functions in time. 

Iterative Solution 

Because the equivalent stiffness matrix k e is a function of the unknown modal 
displacements, the solution takes an iterative form. The time variation of the 
modal displacements and forces may be expressed as: 

«,(») = lie'"-' /,« = £//”■' (8) 

n n 

where ( A ) indicates the dependency on (0 n . Applying |(8)]to (5) and writing in 
iterative form gives: 

q m = [H m ~ x ]f (9) 

where m is the iteration number and 

[/r- 1 ] = [-colli] + ico n [2C,.co r ]+ [co 1 ] + lock';-' + /3k;- 2 ]] 1 (io) 
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The introduction of the weightings a and /? are to aid in the convergence of 
the solution, with the condition that a + J3 - 1 . 

For stochastic excitation, (9) is rewritten as: 

ip , nT 

( 11 ) 


and 


E[q,q,Y = 


The diagonal elements of [s™] are the variances of the modal displacements. 
For the first iteration, [fcj is zero yielding the covariance matrix E[q r q s ~\ of the 
linear system. For subsequent iterations (m> 1), is determined from (6) 


as: 


[K\ = 


dy 

dq 


9/i 


9/l 
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( 12 ) 


where (from|(7)[), 


E 

dy 1 


dq 
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Recall a zero-mean response is assumed, i.e. E[q\ = 0, which reduces |( 13) 

L 


(13) 

to: 


dy 

dq 




j*=j 

L 


'L h )k l E [q ] q l ] 

i,k=j 


j.k=j j’k=j 

The iterations continue until convergence of the equivalent stiffness matrix 
such that 


m — 1 


(14) 




< £ 


The value of £ typically used is 0.1%. Following convergence, the NxN 
covariance matrix of the displacements in physical coordinates is recovered 
from 

£[x,X.] = 0£[ fw ,]0 r (15) 
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and root-mean-square values are the square roots of the diagonal terms in (15) 


Further post-processing to obtain power spectral densities of displacements, 
stresses, strains, etc., may be performed by substituting the converged 
equivalent stiffness matrix into |(5)| and solving in the usual linear fashion. 

Implementation 


The EL procedure as outlined above was recently implemented within the 
context of MSC/NASTRAN using the DMAP programming language [10]. 
The implementation entails first perfonning a normal modes analysis (solution 
103) to obtain the modal matrices, from which a subset of L modes are chosen. 
Since the form of the nonlinear stiffness tenns are not explicitly known, the 
nonlinear stiffness coefficients a' jk and b r jkl must be determined numerically. 
To accomplish this, a series of inverse problems are performed by prescribing 
displacement fields as linear combinations of modes to the linear static 
(solution 101) and nonlinear static (solution 106) solutions. The nonlinear 
stiffness coefficients are then determined from the resulting linear and 
nonlinear nodal forces. The process by which this is done is covered in detail 
in [8] and [9]. The iterative solution is performed within a standalone DMAP 
alter, which has as its output the root-mean-square displacements in physical 
coordinates, the cross covariance in modal coordinates, and the sum of the 
linear and equivalent linear modal stiffness matrices. The latter may by 
substituted for the linear modal stiffness in the modal frequency response 
analysis (solution 1 1 1) for post-processing. 


NUMERICAL SIMULATION APPROACH 


A numerical simulation analysis was performed to generate time history 
results from which response statistics could be calculated. The particular 
method used was finite element based with the integration performed in 
physical coordinates. Two different integration methods were used depending 
on the response level. For lower level responses, an implicit integration 
method was taken which allowed for larger time steps compared with the 
explicit method. Because the implicit scheme is unconditionally stable, a 
convergence study on the time step was undertaken at each response level to 
ensure adequacy of the time step. This was done by halving the At used until 
the time history response over the calculation period was unchanged. As the 
response level became higher at the higher excitation levels, the At required to 
obtain a converged solution became smaller. When the At required was on 
the order of the time step required for the explicit method, it became more 
efficient to perfonn the explicit integration. Both methods produced identical 
results, so the particular choice of implicit versus explicit was dictated solely 
by the time required to run each analysis. In both the implicit and explicit 
methods, the nonlinear defonnation was handled using a corotational scheme 
[11]. The program NONSTAD [12] was used to generate the numerical 
simulation results. 


In Structural Dynamics: Recent Advances, Proceedings of the 7 th International 
Conference, Vol. II, Ed. Ferguson, Wolfe, Ferman, Rizzi, 2000, pp. 833 - 846. 



Loading Time History 

The time history of the load was generated by summing equal amplitude sine 
waves with random phase within a specified bandwidth using a discrete 
inverse Fourier transform. This generated a pseudo-random time history with 
a specified amplitude and period T. The period was specified by 2" At. In 
this study, At = 50 jus and n = 16, giving a period of 1.63845 in duration. 
The selected At corresponded with that needed for the implicit integration 
scheme used for the lower loading levels. The explicit integration scheme 
used for the higher load levels interpolated the signal between points such that 
the specified loading at each At interval was maintained. A radix-2 number of 
time history samples was chosen to facilitate use of radix-2 FFT algorithms 
employed for the subsequent analysis. An ensemble of time histories was 
generated by specifying different seeds to the random number generator. 

A typical time his tory corresponding to an inertial load of 0.05g RMS is 
shown in Figure 1 . The corresponding probability density function is also 
shown with the Gaussian distribution. The power spectral density for 10 
ensemble averages gives a spectrum level of 1.67 x 10 g“/Hz over a 1500 Hz 
bandwidth as shown. A sharp roll-off of the input spectrum practically 
eliminates excitation of the structure outside the frequency range of interest. 

Transient Response Processing 

The structure is assumed to be at rest at the beginning of each loading. An 
initial transient in the structural response is therefore induced before the 
response becomes fully developed. This transient must be eliminated to 
ensure the proper response statistics are recovered. 

The approach taken to eliminate the transient is shown graphically in Figure 2 
Because the loading is pseudo-random, it is possible to apply multiple periods 
of the load and generate the same statistics for periods of length T beginning 
at any point in time. For example, the statistics of the load are the same for 
the period 0 - T as they are for the period T/4 - 5T/4. In a like manner, for 
the linear condition, the response statistics are the same for any period T 
following the initial transient. Therefore, by computing the desired statistic 
for a moving window of period T, it is possible to identify a point in time t 0 
after which the response statistics do not change significantly. In the present 
study, the root-mean- square response was monitored and a time of t a = 0.55 
was found suitable. For each loading history, a total response of 
T + 0.55 = 2.13845 was calculated and the first 0.55 was discarded. Note that a 
similar argument for the nonlinear condition does not hold because the 
nonlinear response is not periodic over T. Nevertheless, the above approach 
was employed with satisfactory results. 
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Response Statistics 

Response statistics were generated from an ensemble of N= 1 0 time histories at 
each load level. Estimates of the displacement root-mean-square served as the 
basis for comparison with the EL method, which essentially had the RMS as 
its basic unknown. Additionally, confidence intervals for the mean value of 
the RMS estimate were generated to quantify the degree of uncertainty in the 
estimate [13] using: 



st , 

n;a/n 

Viv 


< JU X <x + 


St , 

n;a/n 

y JN J’ 


n = N - 1 


where x and s 2 are the sample mean and variance of the RMS estimates from 
N ensembles, and t n is the Student t distribution with n degrees of freedom, 
evaluated at a / 2 . For the 90% confidence intervals calculated, a = 01. 

Estimates of the displacement mean, skewness, and kurtosis were also 
computed to help ascertain the degree to which the assumptions made in the 
development of the equivalent linearization method were followed. Power 
spectral density and probability density functions of the displacement were 
computed for similar purposes. 

RESULTS 


Validation studies were conducted using an 18-in. x 1-in. x 0.09in. (/ x w x /?) 
clamped-clamped aluminum beam with material properties: 

E = 10.6 x \0 6 p si, G = 4.0 x10 V/, p = 2.588 x KT 4 ^^- 


The beam was subjected to an inertial loading over a computational bandwidth 
of 1500 Hz, as shown in Figure 1 This bandwidth was in excess of the 


desired bandwidth of 1000 Hz to allow for the contribution of the higher order 
modes in the EL solution. This is especially important for capturing the anti- 
resonant behavior. The numerical simulation analysis utilized the same 
loading for consistency, although the effect of the higher order modes are 
automatically realized since the computations are perfonned in physical 
coordinates. Since the loading was uniformly distributed, only symmetric 
modes were included in the analysis. In general, any combination of 
symmetric and non-symmetric modes may be included. 


The NASTRAN model used in the EL analysis was comprised of thirty-six !4- 
in. long CBEAM elements. The EL analysis used a four-mode solution 
comprised of the first four symmetric bending modes. Damping was chosen to 
be consistent with the mass-proportional damping of the numerical simulation 
analysis and at a level sufficiently high so that a good comparison could be 
made at the peaks of the PSD. This dictated a critical damping of 2.0%, 
0.37%, 0.15%, and 0.081% for the first four symmetric modes. The finite 
element model used in the numerical simulation analysis was also comprised 
of thirty-six Vi-m. long beam elements. Both EL and simulation finite 
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element models were checked for convergence by running additional analyses 
with models consisting of 14-in. elements. 

Analysis was perfonned at load levels of 0 .8, 3.2, 12 .8, 51.2, and 204.8 g 
RMS, giving a dynamic range of 48 dB. Figure 3 shows the normalized 
RMS out-of-plane (w) deflection at the beam center as a function of load level. 
The numerical simulation results are shown with 90% confidence intervals of 
the RMS estimate. At the lowest load level of 0.8 g RMS, the response is 
linear as can be seen by the comparison with results from a strictly linear 
analysis (NASTRAN solution 111). A small, but noticeable, difference 
between the linear and nonlinear response is noted at the 3.2g RMS load level. 
The degree of nonlinearity increases with load level, as expected. At the 
highest load level, the nonlinear response calculations predict a RMS center 
deflection of 2.5 times the thickness compared with the nearly 11 times the 
thickness from the linear analysis. 

In order to gain greater insight into the nonlinear dynamics, plots of the time 
history, PSD, and PDF ar e shown for three load levels, 0.8, 51.2 and 204. 8g 
RMS, in Figure 4 through Figure 6| , respectively. Data in the time history and 
PDF plots correspond solely to numerical simulation results. Data in the PSD 
plots correspond to numerical simulation and EL results, where the EL results 
were generated by running a linear analysis (NASTRAN solution 111) using 
the equivalent linear stiffness generated by the EL process described above. 
Also shown in the figures are plots of the nonnalized RMS deflection shape 
for both numerical simulation and EL analyses. 

Results for the 0.8g RMS exci tation lev el are shown in Figure 4 . This 
excitation level was shown (see Figure 3 ) to produce a linear response. As 
expected , the PDF mimics the normally distributed PDF of the input shown in 
Figure 1 The averaged PSD and nonnalized deflected shape show excellent 
agreement between the EL and numerical simulation results. This agreement 
helps to establish the confidence in making comparisons between these two 
fundamentally different analyses. 

Figure 5 shows a nonlinear response associated with the 51.2g RMS excitation 
level. The time history of the center displacement has a visibly higher peak 
probability and the PDF exhibits a flattening at the peak. The PSD from both 
the numerical simulation and EL analyses both show the shifting of peaks to 
higher frequencies compared with the linear solution. This shifting is 
associated with a hardening spring type of nonlinearity. Note that the 
frequencies associated with the peaks of the EL solution are the natural 
frequencies of the equivalent linear system. The numerical simulation results 
correctly show the peak broadening effect, which the EL analysis is unable to 
capture. Also distinguishable in the numerical simulation results are 
contributions of hannonics. The effect of nonlinearity is somewhat over 
predicted in the EL result, as seen in the normalized deflection shape. 


In Structural Dynamics: Recent Advances, Proceedings of the 7 th International 
Conference, Vol. II, Ed. Ferguson, Wolfe, Ferman, Rizzi, 2000, pp. 833 - 846. 





The highest degree of nonlinearity is shown in [Figure 6. corresponding to the 
204. 8g RMS load. The time history is further peak oriented and the PDF is 
nearly flat. The peak broadening in the PSD of the numerical simulation 
results is severe, and nearly flattens the spectrum above 350 Hz. Of particular 
interest is flattening of the nonnalized deflection shape, which both the EL 
and numerical simulation results reflect. This is likely to have a significant 
effect on the strain distribution. It is unclear why the results at the 204. 8g 
RMS load level compare more favorably than those at the 51.2g RMS level. 
This behavior warrants further investigation. 

Moments of the center displacement were calculated fro m the n umerical 
simulation results for all load levels. They are provided in Table 1 with the 
RMS center displacement from the EL analysis. The EL and numerical 
simulation results agree well, thus validating the EL analysis over a substantial 
load range. The validity of assumptions made in the development of the EL 
method are ascertained by observing the mean, skewness and kurtosis. The 
mean value is effectively zero for all load levels, indicating the assumption of 
zero mean response has not been violated. Although the PDF is more or less 
skew-symmetric, the shape is flattened at the higher load levels as indicated by 
a decreasing kurtosis from the linear value of 3. The decreasing kurtosis 
values indicate a violation of the Gaussian response assumption. However, 
even with this non-Gaussian response distribution, the EL analysis gives a 
good prediction of the RMS response. 

Table 1: Moments of the center displacement. 


Load Mean (in.) EL Numerical Simulation RMS Skewness Kurtosis 
(§)rms RMS (90% Confidence Interval) 

(in .) (inO 


0.8 -6.54x10 

0.0038 

0.00379 < RMS <0.00381 

0.0197 

3.05 

3.2 -5 .26x1 O' 6 

0.0147 

0.01438 < RMS <0.01532 

0.0190 

2.82 

12.8 -1.62x1 0' 5 

0.0458 

0.04680 < RMS < 0.05350 

-0.0017 

2.47 

51.2 -4.93xl0' 6 

0.1078 

0.11268 < RMS <0.12325 

-0.0002 

2.25 

204.8 3. 78x1 O' 4 

0.2285 

0.21654 < RMS <0.23216 

-0.0058 

2.29 


DISCUSSION 

The RMS random response predictions from the EL implementation have been 
validated through a wide range of load levels. Comparisons with numerical 
simulation results are good, even when the assumption of Gaussian response 
has been violated. 

Differences that do exist warrant some discussion. It is seen that the EL 
approach slightly over-predicts the degree of nonlinearity compared to the 
numerical simulation results. This does not appear to be due to a violation of 
the assumption of a Gaussian response because the over-prediction does not 
correlate with increasing kurtosis of response. The likely reason for the 
difference is in the error minimization approach used. Results previously 
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computed using the alternative strain-energy error minimization technique, not 
included here, indicated better comparison with the exact F-P-K solution than 
the conventional force-based error minimization undertaken here [9]. The 
improved comparison of strain-energy error minimization and simulation 
results was also seen in a different EL implementation for the problem of a 
beam on an elastic foundation [6]. The strain-energy analysis is beyond the 
scope of this paper and is left as an area for further investigation. 

Some implications on the use of the EL technique as a basis for fatigue life 
calculations are worth mentioning. First, assuming that stresses or strains 
from the EL technique will compare equally well with the fully nonlinear 
results, a simple fatigue-life calculation based on RMS levels will be much 
less conservative than calculations based on linear analyses. This offers the 
potential for substantial weight savings for structures designed using nonlinear 
methods. Secondly, it appears that a nonlinear analysis, EL or otherwise, is 
required to accurately calculate the RMS deflected shape. Use of a linear 
RMS deflected shape scaled to the nonlinear level would inaccurately reflect 
the spatial distribution. Simple fatigue-life calculations based on the RMS 
stress or strain could be significantly affected as these quantities depend on the 
spatial distribution of the deformation. Lastly, use of the EL derived PSD 
response in a more sophisticated fatigue-life calculation requires careful 
investigation. Recall that peaks in the equivalent linear PSD may occu r at 
different frequencies than the fully nonlinear PSD, as shown in Figure 5 and 
Figure 6 Methods such as spectral fatigue analysis [14], which take moments 
of the PSD, may incorrectly account for the contribution a particular frequency 
component in the cycle counting scheme. It is not known, for example, if the 
narrowly shaped, higher fundamental frequencies of the equivalent linear PSD 
result in conservative or non-conservative estimates of fatigue life relative to 
predictions made using the fully nonlinear PSD with more broadly shaped, 
lower fundamental frequencies. An assessment of this effect is left as an area 
for further study. 

The question of computational efficiency has not been addressed in this paper 
because the differing analyses were perfonned on different computer 
platforms. What can be stated is that the computational burden for the EL 
approach will increase only slightly with an increase in the size of the physical 
model for the same number of modes used in the solution. This increase is 
associated with the solution of a larger linear eigenvalue problem. The 
expense of the numerical simulation solution, however, will increase 
dramatically with an increase in physical system size. 
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Figure 1: Typical loading time 
history, probability density, and 
power spectral density. 



Time (s) 

Figure 2: Transient response 
processing. 



Figure 3: Normalized RMS center 
deflection as a function of inertial 
load 
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Figure 4: Center deflection 
response and normalized deflection 
shape at 0.8g RMS inertial load. 


Figure 5: Center deflection 
response and normalized deflection 
shape at 5 1 ,2g RMS inertial load. 
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Figure 6: Time history, PSD, and PDF of center deflection response and 
normalized deflection shape at 204. 8g RMS inertial load. 
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